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Abstract 

We derive the generalised master equation for reaction-diffusion on networks from an underlying 
stochastic process, the continuous time random walk (CTRW). Using this model we have investi- 
gated different types of pattern formation across the vertices on a range of networks. Importantly, 
the CTRW defines the Laplacian operator on the network in a non ad-hoc manner and the pattern 
formation depends on the structure of this Laplacian. Here we focus attention on CTRWs with 
exponential waiting times for two cases; one in which where the rate parameter is constant for all 
vertices and the other where the rate parameter is proportional to the vertex degree. This results 
in non-symmetric and symmetric CTRW Laplacians respectively. In the case of symmetric Lapla- 
cians, pattern formation follows from the Turing instability. However in non-symmetric Laplacians, 
pattern formation may be possible with or without a Turing instability. 
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I. INTRODUCTION 



Networks have been extensively studied as models for highly connected systems in biology 
[1] , physics [2] and the social sciences [3] . Here we consider the dynamics of reaction-diffusion 
networks with reactions of particles on vertices and diffusion of particles between vertices. 
One of the interesting features of these systems is the possibility of spontaneous pattern 
formation whereby the concentration of particles is non-homogeneously distributed across 
the network. 

React ion- diffusion models form the basis of pattern formation in a wide range of self- 
organising systems. Turing's seminal work in 1952 jl] showed that in a two species system 
with reaction and diffusion, patterns may arise (Turing patterns) due to an instability in 
the reaction dynamics caused by differing rates of diffusion. Such patterns emerge in bio- 
logical morphogenesis [SHI] chemical reactions [lOj, propagation of viruses [HI [12] as well as 
ecosystems encompassing competing animals |13j . 

In early work, the existence of Turing patterns on networks was explored by Othmer and 
Scriven [H] whereby they showed that the structure of the network has a direct effect on the 
resulting pattern. Similar behaviour was also reported by Colizza et. al. Subsequent 
studies has considered scale-free networks \16\, coupled reactors [17] and functional gene 
networks [18]. 

Recent studies of reaction and diffusion on networks have identified that multiple coex- 
isting stationary states exist for many systems [T9H22] . The effects of feedback [20] and the 
formation of travelling fronts [22] have also been studied in these systems. 

A discrete space description of diffusion can be modelled by a Random Walk (RW) [231 - 
[2S]- RW on networks have been extensively studied in this context. See for example [2SH2H]- 

The continuous time random walk (CTRW) is a generalization of the random walk in 
which the random walker waits a time (drawn from a waiting time probability density) before 
jumping [291 ED] • This model has been particularly useful to model diffusion in systems with 
disorder in waiting times resulting in anomalous diffusion |31j- The continuous time random 
walk has been further generalized to include linear reactions [321 133] , linear reactions with 
multiple species [33] and nonlinear reactions [35H39] . 

Continuous time random walks on networks without reactions have been studied in [IQ] . 
Here we consider continuous time random walks with reactions on networks. The descrip- 
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tion is at the mesoscopic level, with particles jumping between vertices and with reactions 
between particles at vertices. This model description could be applied to the analysis of dif- 
fusion tensor imaging data jHJ |32] and to spatio-temporal models in epidemiology [151 HI] . 
The essential assumption in this description is that there is an underlying stochastic process 
that defines the motion of particles on the network that may be modelled as a CTRW. 
Specifically, we assume that each vertex in the network can be occupied by any number of 
walkers and the reactions among the walkers at a vertex are governed by the same reaction 
kinetics at all vertices. Our principal concern is then the collective motion of the walkers on 
the network as a whole. 

Within this CTRW framework we have derived a family of diffusive network Laplacian 
operators. The reaction-diffusion behaviour with these network Laplacian operators may 
differ from that with the continuum Laplacian operator In the following we have 

considered pattern formation arising from continuous time random walks on networks with 
reactions using the different forms of the Laplacian. 

The remainder of the paper is as follows. In section [TT] we derive the master equations 
that describe the continuous time random walk network reaction diffusion model. In section 



III we describe different pattern formation mechanisms that may arise depending on the 



form of the CTRW Laplacian. In section IV we present numerical simulations of pattern 
formation in the network models. We conclude with a summary and discussion in section [V] 



II. DERIVATION OF A CTRW NETWORK REACTION DIFFUSION MASTER 
EQUATIONS 

Diffusion, or diffusive-like phenomena, can arise on a network from a variety of source 
depending on the specifics of the network. In considering diffusion we should begin by 
defining a stochastic process that make physical sense for the phenomena on the network. 

A continuous time random walk is a stochastic process that naturally limits to diffusion 
in the continuum [211 EQj- To model a reaction-diffusion process we assume that the motion 
of each individual particle can be expressed as a CTRW. That is to say the particles will 
jump from vertex to vertex on the network according to the edges present. On each vertex 
they will wait for a random time before randomly jumping to a connected vertex. This 
is analogous to spatially varying diffusion that has been previously studied in the spatial 
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continuum case |16] . 

The reactions occur between particles that occupy the same vertex. We consider the 
reactions to be a birth-death process of the particles. Particles will be created according to 
some probability, and destroyed according to a different probability. These probabilities may 
depend on the density of other particles on the vertex. We can then derive the equations 
that govern the evolution of a single particle in a time. In this manner the evolution of a 
single particle is subject to only the death probabilities. The birth process is included by 
summing the initial conditions of each single particle CTRW with the probability that a 
particle was created on a particular vertex at an instant in time. 

The assumption that each vertex is a well mixed system is made so that the probability 
of a particle being destroyed in reactions is not dependent on the amount of time a particle 
has been waiting on the vertex. In a similar manner we also assume that each newly created 
particles waiting time is independent of the parent particles waiting time, similar to Model 
B in [38j which was first considered by Vlad and Ross [35]. The first step in this process is 
the derivation of the master equation for the evolution of a single particle subject to a time 
and space inhomogenous probability of death. 

A. Single Particle CTRW Death Process Density 

Consider a graph whose vertices form the set W = {wi, wj} where J is the number of 
vertices. Let p{wj,t\wo,0) be the probability density for a random walker to be on vertex 
Wj at time t given it started on vertex wq E W a.t time t = 0. 

Define qn{wj,t\wQ,0) as the conditional probability density for arriving at vertex wj at 



that a particle stays alive from t' to t given it does not leave vertex Wj where /3 is a death 
rate that in general may depend on vertex and time. The initial condition for n = is given 



f /S{wj,t")dt' 



time t after n steps. We define the reaction survival function, e 



as the probability 



by 



(1) 



In general, we can write 




(2) 
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where '^{wi,Wj,t,t') is the probabihty density of the transition to vertex wj at time t given 
the random walker arrived at vertex Wi at an earher time t' after n steps. 

We assume that '^{wi^Wj,t,t') may be defined as a product of the independent jump 
probabihty density \{wi,Wj) and tpiwjjt). 

"^{wi, Wj,t, f) = X{wi, Wj)ip{wj,t - t'). (3) 

The Wj in ip acknowledges that the waiting time density may be a function of the vertex. 
The conditional density for the walker to arrive at wj at time t after any number of steps is 
found by summing over all n steps using Eq. Q and Eq. g. 

oo 

q{wj, t\wo, 0) = ^ qn{wj, t\wo, 0) (4) 

n=0 

°° * f I3{wj,t")dt" 



n=0 i=\ n 



'{wi,Wj,t,t')e qn{wi,t'\wo,0)dt' 

(5) 

/■ - f I3{wj,t")dt" 



1=1 



(6) 

We can then define the conditional probability density for the random walker to be at 
vertex Wj at time t; 

t t 

r - J i3{wj,t")dt" 

p{wj,t\wo,0) = I (j){wj,t — t')e *' q{wj,t'\wQ,Q)dt' , (7) 



where 4>{wj,t — t') is the probability that the particle does not jump during the period of 
time t — t'; 

t-t' 

(piwj, t-t') = l- J ip{wj, t")dt". (8) 



B. Single Particle CTRW Death Master Equation 

The derivation of the master equations describing a CTRW death process on a network 
is similar to the derivations presented in [17], [18] and [38] • Formally, the integrals over 
probability densities should be treated as Riemann-Stieltjes integrals and care has to be 



taken due to the discontinuity in the arrival density q{wj,t\wo,0) at time t = [38]. To do 
this, write 

q{wj, t\wo, 0) = 5^^^^^S{t - 0+) + q^{wj,t\wo, 0) (9) 

where g"*" is right side continuous at t = 0. Thus by substitution of Eq. ^ into Eq. ([t]) we 

get 

t t t 

p{wj,t\wo,0) = Suj^^wo4>{wj,t)e + / q^{wj,t'\wo,0)e *' (j){wj,t - t')dt' . 



(10) 

We now differentiate this equation with respect to time, using the Leibniz rule for differ- 
entiating under the integral sign, to obtain 

t 



dp(wi,t\wo,0) , X /■ /, -Ifii^j,t")dt" 

' ' = g+K-,tko,0) - / g+K-,tVo,0)e ijiwj,t - t')dt' 



- J I3{wj,t")dt" 

/3{wj, t)p{wj,t\wo, 0) - 5^„^,^oe o ipiwj, t). 



(11) 



Define the flux leaving vertex Wj at time t as 

t t 
- f p{wj,t')dt' r „ „ , 

i{wj,t\wo,0) = ^uj^^wo^ ° tlj{wj,t) + / q^{wj,t'\wo,0)e tlj{wj,t - t')dt' 





- f I3{wj,t")dt" 



- J I3{wj,t")dt' 

q{wj,t'\wo,0)e y^^-j 



ip{wj, t — t')dt' . 



(12) 
(13) 



We can then rewrite Eq. (11) as 

dp{wj,t\wo,Q) 



dt 



g+(wj,t|wo,0) - i{wj,t\wQ,Gi) - P{wj,t)p{wj,t\wo,0). (14) 



Using Eq. ([6]), ^ and (13), the rate of arrivals at vertex Wj can be expressed as 

J 

q^{wj,t\wo,0) = \{wi, Wj)i{wi, t\wo, 0). 



(15) 



i=l 



Now we can define p through an evolution law as follows 



dp{wj,t\wo,0) 
dt 



X{wi,Wj)i{wi,t\wo,0) - i{wj,t\wo,0) - (3{wj,t)p{wj,t\wo,0). (16) 
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Following Fedotov [38] we can find an expression for the flux i in terms of p using Laplace 
transform methods on Eq. ([T]) and Eq. ( 13 ) respectively. We first divide both equations by 

- J I3{wj,t")dt" 



this yields; 



C < p{wj, t\wo, 0)eo 



J I3{wj,t")dt" 



J I3{wj,t')dt' , 

C{q{w„t\wo,0)eo } C{(j){wj,t)} (17) 



and 



L < i{wj, t\wo, 0)eo 



f li{Wj,t")dt" 



f I3{wi,t')dt' 

C(q{wj,t\wo,0)eo } C{ij{wj,t)}. (18) 



Rearranging Eqs. (17) and (18); 



C < i{wj, t\wo, 0)eo 



J I3{wj,t')dt' 



C < p{wj, t\wo, 0)eo 



Inverting the Laplace transform, we get; 



i{wj,t\wQ,0) = / K{wj,t — t')p{wj,t'\wo,0)e *' 



J I3{wj,t")dt'' 



dt' 



where the memory kernel is defined by; 



^ ^' ' \C{<j>{w,,t)}j 



(19) 



(20) 



(21) 



The master equation for the CTRW process on a network is found by the substituting Eq. 



(20) into Eq. (16) 



dp{wj,t\wo,0) 



dt 



y^X{wi,Wj) / K{wi,t ~ t')p{wi,t'\wo,0)e *' dt'- 



- f I3(wj,t")dt" 

K{wj,t — t')p{wj,t'\wo,0)e *' dt' — (3{wj,t)p{wj,t\wo,to). 



(22) 



In the case of power law waiting time densities on a uniform grid network, this is similar to 
Eq. (29) in [39j. 

C. Ensemble CTRW Birth-Death Master Equations 



To describe a birth-death process we also need to account for the creation of new particles, 
and hence need to consider ensembles of particles. We define ri{wj,t) as the probability of 
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a particle being created at vertex wj and at time t. Then we can define 



uiw 



J' 



t 

t) = j p{wj,t\wo,to)r]{wo,to)dto 



(23) 



as the density of a particles at vertex Wi at time t. 



Taking care to differentiate Eq. (23) using Leibniz rule, we then substitute in Eq. (22) 
and simplify to get 

t 

du{wj, t) 



dt 



I \ I I \ / / N f^P('W^7; ^I'W^O; ^o) , 

ri[wQ,t)p[Wj,t\wQ,t)^ \ r]{wo,to) — dto 



E 

WO&W . 



dt 



(24) 



1=1 



fP{w,,t")dt" 



/ K{wi,t -t')X{wi,Wj)e *' p{wi,t'\wo,to)dt'- 



- f I3{wj,t")dt" 

K{wj,t — t')p{wj,t'\wo,0)e *' dt' — (3{wj,t)p{wj,t'\wo,tQ) 



dto + v{wj,t). 
(25) 



As p{w, t\wQ, to) = for all t < to we can write; 



du{wj, t) 
dt 



t 


" J 


I 







.1=1 


- f I3{wj,t")dt" 

e 



■ f I3{wa")dt" 



^ / p(wi,tVo,^o)'7(w^o,^o)c?^o- 



K{wj,t-t') ^ / p{wj,t'\wo,to)T]{wo,to)dtQ 



dt'- 



l3{wj,t) ^ / p{wj,t\wo,to)r]{wo,to)dto + r]{wj,t) 



Finally we arrive at the master equation for a CTRW with reactions on networks 



(26) 



du{wj, t) 
dt 



t r 



>- 



-jp{wa")dt" 

y^K{wi,t — t)X{wi,Wj)e *' u{wi,t)- 



1=1 



(27) 



- J I3{wj,t")dt" 

K{wj,t-t')e *' u{wj,t') 



dt' — f3{wj,t)u{wj, t) + ri{wj, t). 



which can be written in the general form 

du{wj, t) 



dt 



L[u{wj,t)] + f{u{wj,t)) 



(28) 



where 



L[u{wj,t)] 



K{wi,t — t')X{wi,Wj)e *' u{wi,t') 

.1=1 



- f l3(wj,t")dt" 

K{wj, t — t')e *' u{wj, t') 



dt' 



is the CTRW network Laplacian and 



f{u{wj,t)) = -(3{wj,t)u{wj,t) + r]{wj,t). 



(29) 



(30) 



D. CTRW Laplacian with Exponential Waiting Times 



We will now apply Eq. (27) to the case of exponential waiting times, 



^{wj,t) = a(«;j)e-°("'^)*. 



(31) 



This greatly simplifies the master equation. In general, the Laplace transform of the waiting 
time density, ipiwj.s) = t)} and the Laplace transform of the survival probability 

(f){wj, s) = C{(j){wj,t)} are related by (piwj, s) = 



As ipiwj^s) is exponential, then ip{wj,s) = so 
Laplace transform a{wj)6{t). Thus 

K{wj,t) = a{wj)6(t). 



(f>(Wj,s) 



a(w^) with the inverse 



(32) 



Thus by substitution, we can rewrite the master equations Eq. (27) as 



= ^ a{wi)X{wi, Wj)u{wi, t) - a{wj)u{wj, t) + f{u{wj,t)) 

1=1 
J 

= ^ Liju{wi, t) + f{u{wj, t)) 
1=1 

where L is a member of the family of general CTRW network Laplacians defined as 

Lij = {a{wi)X{wi,Wj) - a{wj)6ij} 



(33) 
(34) 



(35) 



In the following we consider two special cases of Eq. (35 ). For both cases, we assume that 



the jump probability and the waiting time probability are only functions of vertex degree 
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and time respectively. Both of these cases have been previously considered without reactions 
To describe the network we used the adjacency matrix 

, 1 if vertices Wi and Wj are connected , , 

A,j = { ' (36) 

otherwise. 

Each vertex, Wj, has a degree, kj, that is defined as the total number of edges that link to 
the vertex. This is also the sum of each row of the adjacency matrix. 



1. Case A 



We assume that the waiting time on each vertex is identical and that the probability of 
jumping across an edge is equal for a given vertex. Formally, we let a{wj) = for all 
Wj G W and let \{wi,Wj) = ^. Thus our Laplacian becomes 

= aA Qr.^iJ - kj^ (37) 

Note that in this case the Laplacian is not symmetric. The steady state for a CTRW with 
no reactions using this Laplacian, in general, will not be uniform on the vertex set. 



2. Case B 



An alternative is to let the waiting time on each vertex change proportionally to the 
vertex degree. This allows the rate of particles jumping along each edge to be constant. 
Formally, we let a{wj) = askj for all Wj G W and let X{wi,Wj) = ^. Thus our Laplacian 
can be described as 

Li J = as (Ajj - kj5ij) (38) 

This is the well studied graph Laplacian [IS] [20] • It is important to note that if a connected 
network is regular, i..e kj is the same for all Wj, the cases are equivalent up to a scale factor. 



III. PATTERN FORMATION 



In continuum reaction-diffusion partial differential equations the concentration may vary 
on the spatial domain. This patterning also holds on discretisation of the spacial manifold 
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where the Laplacian operator is replaced by a discrete Laplace Beltrami operator [50] [5T] . 
This can be considered as a special case of pattern formation on a discrete network. In 
general variation in concentrations on each vertex may occur on any network. This variation 
may arise in a number of different ways. Firstly if the Laplacian matrix is not symmetric in 
a system without reactions, there will be a build up of concentrations on vertices according 
to their degree. This steady state pattern will be a multiple of the eigenvector of the 
Laplacian with zero eigenvalue. Secondly with a non-symmetric Laplacian matrix in systems 
where the reactions have a finite non-zero steady state solution, the interplay between the 
diffusion and the reactions will cause a different pattern across the network. Unlike the no 
reaction case, in this pattern vertices with the same degree may have different concentrations. 
We refer to these two mechanisms as Laplacian pattern formation, as the patterning is 
driven by the non-symmetric Laplacian matrix. Lastly if the reaction terms permit a Turing 
instability, whereby the spatially homogeneous steady state becomes unstable in the presence 
of diffusion, then a Turing Pattern may form |4j. This instability is permitted for both a 
symmetric and non-symmetric Laplacian matrix. 



A. Laplacian Pattern Formation 

To completely eliminate any interplay with Turing patterns, we will consider a single 
species model that cannot permit a Turing instability [H2]. 

^ = L.u + f (u) (39) 

where u = u{wi,t), u{wj, t) ) and f (u) = f{u{wut)), f{u{wj, t)) ) • 

If we take the trivial reaction term, f(u) = 0, the only possible form of spatial pattern 
formation comes from a non- symmetric Laplacian. The master equation is then simply: 

f^..u (40) 

It is clear that if L is non-symmetric, there will be a nonuniform rate of particle transport 
along the network that produces a pattern of different concentrations of particles in the 
long term. This concentration vector must be a multiple of the eigenvector of the Laplacian 
matrix with eigenvalue zero as it is a solution of the the vector u that satisfies the steady 



state of Eq. (|40j), i.e. 

L.u = 0. (41) 
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If the reaction term f (u) 7^ and has a finite non-zero equihbrium solution then the pattern 
will no longer correspond to an eigenvector of the Laplacian matrix. 

Define the reaction steady state u* such that the reaction term f (u*) = 0. By considering 



the behaviour of u = u* + Au where Au is some perturbation, Eq. (39) may be linearised 
to become 

^ = Df(u*).Au + L.u* + L.Au. (42) 

The steady state solutions are given by 

Au = -{L + Df{u*)y\L.u*, (43) 

provided the inverse exists where Df{u*) is the derivative operator of the reaction vector 
field f at u*. This is the linear prediction of the pattern. In the case where the reaction 
term is linear in u this equation is the exact solution for the pattern. It should also be 
noted that in the case of a symmetric Laplacian L.u* = 0, as the constant vector is always 
a null vector of a symmetric Laplacian. This equation does not hold without reactions as a 
non-symmetric Laplacian is in general non-invertible. 

This type of pattern formation is very different from a Turing pattern. A Turing pattern 
is formed by an instability in the dynamics whereby an otherwise stable solution is made 
unstable by the presence of diffusion. It is a bifurcation phenomena as the diffusion must 
reach a critical value before the solution becomes unstable. In this case however we have 
a pattern that will form for any amount of diffusion, the stability of the solution does not 
change but rather the solution itself is a function of the diffusion. To examine this we 
introduce a scale parameter s for our Laplacian whereby the linear predictor for the pattern 
becomes 

Au = - {sL + Di{u*)y^ sL.u* (44) 

The new scalar parameter that governs the speed of diffusion in the system. In the limit 
s — ^ we have; 

Au = (45) 

which corresponds to no pattern, and each vertex taking the equilibrium values of the 
reactions. In the limit s — > 00 care has to be taken with the existence of the inverse. Taking 



Eq. (44) and multiplying through by {sL + Di(u*)) dividing by s and taking the limit we 
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have 

lim (l + Au = -L.u* (46) 

s-s>oo \^ -5 / 

L.Au + L.u* = (47) 
L.u = (48) 

and so the pattern will be equivalent to the case with no reactions. In this way it can be 
seen that this type of pattern formation is the mixing of the pattern formation from the 
non-symmetric Laplacian and the reaction equilibrium. 

B. Turing Pattern Formation 

Turing pattern formation can arise from a Turing instability Such instabilities were 
originally proposed by Turing as a possible explanation for morphogenisis. We consider 
a general two species reaction-diffusion system with particle concentrations u and v. The 
master equations, [see Eq. (pi))], can be written as 

^'^^'^'^^ = J2auLiju{wi,t) + f{u{wj,t),v{wj,t)) (49) 
1=1 

and 

dv(wj,t) 

Q-^ = 2_^ayLijv{wi,t) + g{u{wj,t),v{wj,t)) (50) 

1=1 

where the functions f{u{wj,t),v{wj, t)) and g{u{wj, t),v{wj, t)) incorporate the creation and 
destruction probabilities, ie; 

f{u{wj,t),v{wj,t)) =r]^{wj,t) - (3u{wj,t)u{wj,t) (51) 
= r]uiu{wj,t),v{wj,t)) - (3uiu{wj,t),v{wj,t))u{wj,t) (52) 

and similarly for g{u{wj,t),v{wj,t)). Here note that the a„ and is factored out of the 
Laplacian, Lij so that the operator is the same in both equations. 

1. Linear Stability Analysis 



To consider Turing instabilities, we first rewrite Eqs. (49) and (50) in vector form 

^ = AX + F(X) (53) 
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where 



X 



u{wi, t), . . . , u{wj, t), v{wi, t), . . . , v{wj, t) ) 
^ Xi, . . . , Xj, Xj+i, . . . , X2J j , 

( f{X^, Xj+,), f{Xj, X2J, g{X,, XJ+^), g{Xj, X^j) 
[f,,...,Fj, Fj 



' J+l! • • • ! F2J 



and A 



auL 
a^L 



Linearising about the steady state X* with F(X*) = and X = X* + AX. We then 



substitute this into Eq. ( 53 ) to get 

aAx 



dt 



(A + DF(X*)).AX + A.X 



where 



DF{X*) 



dF 



dX, 



(54) 



(55) 



We now apply the affine transform AY = AX+ (A + DF{X.*)) ^ .A.X* to Eq. yielding 

dAY 



dt 



{A + DF{X.*)) .AY 



(56) 



which gives solutions of the form 



AY = J2e''''Pj{t) 



(57) 



where Hj is the j eigenvalue of (A + DF(X*)) with corresponding eigenvector uj. Pj 
is a polynomial in t for repeated eigenvalues. The long time behaviour of AY is then 
approximated by 



AY ~ e^*'P*{t)v* 

where /i* corresponds to the real component of the largest eigenvalue. 

In the linear stability analysis the concentrations of the two species evolve as 

X ~ e^**P(t)z/* + X* - (A + DF{1C))~^ A.X* 



(58) 



(59) 
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2. Dispersion Relation 



In the continuum case the real components of the eigenvalues (i.e.: stability) of the 
homogeneous steady state can be plotted against spatial frequency to obtain a dispersion 
relation showing the range of spatial frequencies that will grow with time. In an analogous 
manner for the network case, we plot the stability of the reaction-diffusion system as a 
function of a scale parameter s equivalent to scaling the waiting time for both species on 
the network: 

— = sA.X + F(X) (60) 

When s is zero there is no coupling between vertices and each vertex will be in equilibrium 
according to the reaction equations. We identify two critical values of s greater than zero 
from our linear stability analysis. Firstly the Laplacian type pattern arising from the cou- 
pling of reactions to diffusion may change to a Turing pattern at a critical value of s. As s 
is further increased the Turing pattern, if it occurs, will persist till the second critical value 
when the pattern reverts back to a Laplacian type pattern. 



IV. EXAMPLES OF PATTERN FORMATION 

To illustrate the properties of both the Laplacian and Turing patterns on networks the 
master equations with exponential waiting times were solved numerically. For the Laplacian 



patterns, Eq. (34) was solved with Logistic reaction kinetics using the Case A Laplacian 



operator. The Turing patterns were examined by solving Eqs. (49) and (50) with Gierer- 
Meinhardt reaction kinetics and both the Case A and Case B Laplacian operators. In 
both cases the equations were solved on random networks generated by the Barabasi-Albert 
algorithm [3], and the Watts-Strogatz algorithm [5^. The equations were also solved on 
deterministic binary tree networks. 

The first reaction kinetics that we consider is the logistic equation |54j . 

f {u{x,t)) = ru{x,t){l — u{x,t)) (61) 

where r is a constant. This governs the growth rate of a single species. In the following 
examples we will take r = 1. When applied to a network, this simple example could 
be considered a model for animal populations in a set of connected habitats, where the 
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population is constrained by natural limits. The reaction-diffusion master equation in this 



case is found by substituting Eq. (61) into Eq. (39) 



du(wj,t) sr~^ r I \ I \/i / w 
= 2_^Liju{Wi,t) +u{Wi,t){l -u{wi,t)) 

i=l 



(62) 



Here L is the Case A Laplacian with a = 1, i.e. Lij = — Sij. 

The second model we consider has Gierer-Meinhardt reaction kinetics [3]. This is a two 
species model that permits Turing instabilities. The reaction terms in the model are: 



and 



u(x 

f{u{x, t),v{x, t)) = cp ' - /i u{x, t) + Pop 

f (x, t) 



g{u{x, t),v{x, t)) = CdP u{x, ty — V t>(x, t) 



(63) 



(64) 



We can apply this model to a network setting by substituting Eqs. (63) and (64) in Eqs. 



(49) and (50): 



du{wj,t) u{wj,t) 

~~dt """^ 

dv{wj,t) 



p u{wj, t) + Pop + a« ^ Liju{wi, t) 



i=l 



dt 



In the following we use; 



CdP u{wj, ty — v v{wj, t) + Lijv(wi, t). 



(65) 
(66) 



i=l 



Po = 1, P = 1, ^ = — , P = , c = 1, and Cd = . 

, y , 32' ^ 256' ' 128 



(67) 



For the Case A Laplacian we will use: 



a,, = 1 and a, 



1 

256' 



(68) 



To place the Case B Laplacian, i.e. Lij = Aij — ki6ij, on a comparable footing to the Case 
A we rescale the parameters and Ot, to ensure that the mean waiting time across the 
network is comparable in both cases. Explicitly 



J 



and ay 



1 



J 



256; EU^j 



(69) 
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A. Barabasi- Albert Network 

The Barabasi-Albert (BA) network is a random network with a power law distribution 
of vertex degrees [3]. The network is iteratively generated by adding a vertex at each step 
that connects to k existing vertices where the probabihty of attachment is proportional to 
the degree of the existing vertices. 

We first consider a purely diffusive process governed by the Case A Laplacian on BA 
networks. In this case, the concentration on each vertex is proportional to its degree, re- 
sulting in the shape of the concentrations to be similar to a power law shape as shown in 
Fig. [T| This distribution, which corresponds to the eigenvector of the Laplacian with a zero 
eigenvalue, is monotonic with increasing vertex degree. 

04- 
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u 
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0.1- ^1 - 

o.o k , , , , ,- 

100 200 300 400 500 
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FIG. 1. (Colour Online)Case A Laplacian diffusion on a BA network with A; = 3. A characteristic 
pattern of a network with 50 vertices {left) and the distribution of concentrations on a 500 vertex 
network {right). The vertical lines demark regions of different vertex degree 



O « 

o o 

D OO 

O 



o ^ 
o 



9) 

O Q 



O 
O 

o 

o 
o 



o 



oO o o 



With the addition of the logistic reaction term, concentrations across the network change. 
The overall power law shape is still present, however the pattern is no longer monotonic with 



respect to vertex degree. The linear predictor of the pattern across the network, Eq. (43), 
is a good approximation near the reaction steady state, see Fig. [2j The mean of the 
concentration across the network is reduced from the reaction equilibrium value. 

When the reactions are governed by the Gierer-Meinhardt model, Turing instabilities can 
arise producing patterns different to the Laplacian patterns. We first consider the patterns 
when the Case A Laplacian is used, a representative example is shown in Fig. |3| The 
exact pattern is determined by the initial conditions of the system. This multi-stability is in 
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FIG. 2. (Colour Online) Case A Laplacian with logistic reaction kinetics on a BA network with 
k = 3. Concentration {left) and difference of linear predictor and concentration (right). The 
pale horizontal line is the reaction steady state value and that the data is ordered according to 
concentration within each segment of constant vertex degree. The vertical lines demark regions of 
different vertex degree 



contrast to the Laplacian patterns that have no initial condition dependence. In all observed 
Turing patterns the concentrations for both types of particles are split on vertices with low 
degree. There is also an increasing concentration as a function of vertex degree as seen in 
the previous BA patterns, especially for the concentration of u particles. 

When considering the Case B Laplacian, as shown in Fig. [4], the patterns shown are 
clearly very different to those in Fig. [3j There is a splitting of concentration across all 
vertex degree segments and the lower level shows practically no increase in concentration as 
a function of vertex degree. Moreover the concentrations achieved are higher and lower for 
u and V respectively when compared to Case A. Unlike for the Case A model, the patterns 
exhibited by u and v have similar profiles. 

It is interesting to note the change in mean concentration for the two Turing patterns. 
For the Case A pattern we have both the the mean of u and v greater then the reaction 
steady state value with no diffusion. For Case B, the opposite is true and the mean values 
are less then the reaction steady state values. 

The dispersion relations for Case A and B are plotted in Fig. |5j The Turing instability 
occurs over a larger range in Case B than Case A. 
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FIG. 3. (Colour Online) Case A Laplacian Gierer-Meinhardt reaction diffusion on a BA network 
with k = 3. Concentration of u (top) and v (bottom) on 50 {left) and 500 {right) vertex networks. 
The horizontal line gives the reaction steady state. 



B. Watts-Strogatz Network 

The Watts-Strogatz (WS) network is a random network characterised by the small world 
property of having a low graph diameter |53] . The network is generated by creating a ring 
lattice where each vertex is connected to k adjacent vertices. Then each subsequent edge 
may be reconnected with some probability p to some other vertex chosen uniformly from all 
vertices. 

Once again we first consider a purely diffusive process governed by the Case A Laplacian. 
The concentration on each vertex is proportional to the degree of the vertex as shown in 
Fig. § 

The same network with logistic reaction kinetics is shown in Fig. [7j Similarly to the BA 
network, each segment of vertices with identical degree shows a gentle slope in concentration 
with tapered boundaries. However the underlying concentration distribution is similar to 
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FIG. 4. (Colour Online) Case B Laplacian Gierer-Meinhardt reaction diffusion on a BA network 
with = 3. Concentration of u {top) and v (bottom) on 50 {left) and 500 {right) vertex networks. 
The horizontal line gives the reaction steady state 



that of the previous case. We see that the linear predictor is a good approximation when 
the concentration is near the reaction steady state value. 

We consider the Gierer-Meinhardt reaction kinetics with both Case A and Case B Lapla- 
cian operators as shown in Figs. |8] and [9] respectively. In both cases, the distribution of 
concentrations is bimodal for each vertex degree. 

The example using the Case B Laplacian exhibits very similar patterns to those using 
the Case A. Similarly to the BA network, the highest concentrations in the second example 
are always increasing in contrast to the the first example. Unlike the patterns on the BA 
network, the concentrations for both cases are at comparable quantities. 



Dispersion relations for the Case A and Case B models are shown in Fig. [TO] Unlike in 
the BA network, both Laplacian operators permit dispersion relations with almost indistin- 
guishable profiles. 
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FIG. 5. (Colour Online) Dispersion Relation for Case A (left) and Case B (right) Laplacian 
Gierer-Meinhardt reaction diffusion on a BA network with k = 3. 
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FIG. 6. (Colour Online)Case A Laplacian diffusion on a WS network with A; = 3 and p = 0.1. A 
characteristic pattern of a network with 50 vertices (left) and the distribution of concentrations on 
a 500 vertex network (right). The vertical lines demark regions of different vertex degree. 



C. Binary Tree Netv^ork 

A binary tree is a deterministic network whereby starting with a single vertex, we connect 
two new vertices which are in turn each connected to another two new vertices and this 
process is repeated until the desired number of vertices is reached. 

The concentrations for a system with no reactions and a Case A Laplacian are shown in 



Fig. 11 



There are three distinct concentrations across the network as the concentrations are 
proportional to vertex degree. 

21 




FIG. 7. (Colour Online) Case A Laplacian with logistic reaction kinetics on a WS network with 
k = 3 and p = 0.1. Concentration (left) and difference of linear predictor and concentration (right). 
The pale horizontal line is the reaction steady state value and that the data is ordered according 
to concentration within each segment of constant vertex degree. The vertical lines demark regions 
of different vertex degree. 



With the addition of the logistic reaction kinetics, concentrations are no longer simply 



proportional to vertex degree as shown in Fig. 12 



Steady state concentrations for Gierer-Meinhardt reaction kinetics and both Case A and 



Case B Laplacian operators are shown in Figs. 13 and 14 respectively. 



The dispersion relation of the Case A and B models as shown in Fig. 15 differ significantly 
to those on the BA and WS networks [see Figs. [S] and 10 . 



V. SUMMARY AND DISCUSSION 



In this paper we have derived the generalized master equation for reaction diffusion pro- 
cesses on networks based on the CTRW as the underlying stochastic process. We assumed 
that the vertices of the network can be occupied by many particles with the reactions occur- 
ring amongst the particles on the same vertex. The reactions were assumed to be governed 
by the same reaction kinetics on each vertex. The CTRW models particles jumping between 
vertices and reaction kinetics were incorporated into this as birth and death processes. We 
used the CTRW framework to derive a family of Laplacian operators that govern the dif- 
fusion on the network. These operators are dependent on the waiting time probability 
density at each vertex and the jumping probabihty density between vertices. In the case 
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FIG. 8. (Colour Online) Case A Laplacian Gierer-Meinhardt reaction diffusion on a WS network 
with = 3 and p = 0.1. Concentration of u {top) and v (bottom) on 50 (left) and 500 (right) 
vertex networks. The horizontal line gives the reaction steady state. 



of non-exponential waiting time densities these operators convolve the reaction and trans- 
port processes. However in the case of exponential waiting times the operators are purely 
transport operators. 

In general it is to be expected that the CTRW network reaction-diffusion model can 
lead to unequal concentrations of particles across the vertices. We investigated this pattern 
formation in the concentrations of particles. We considered the jumping probability to be 
equal across any edge from a given vertex and we considered the waiting time densities to 
be exponential, for two choices of the rate parameter; Case A where the rate parameter was 
taken to be proportional to the vertex degree and Case B where the rate parameter was 
taken to be the same for all edges. 

We identified three distinct pattern formation mechanisms in our CTRW network 
reaction-diffusion models. In the case of symmetric network Laplacian operators pattern 

23 



o 



cP°»8 



o °o% 

o oo , 



'OO' 



o 
o 



D 



O' 



O %i 

O OO ^ 



o 



'OO 



o 
o 



I 



I 



.976. 1000 

goo 

U 

1^^'-' 400 


/ 

1 


/ 

I 


I 

i 
■ 


- 30315 200 






^^^^ 












30000 



V 



20000 



200 300 

Vertex 



400 



100 200 300 400 

Vertex 



500 









/ 


A 


1 


/ 

a 

■• / 




f 


^ 1 


t 



500 



FIG. 9. (Colour Online) Case B Laplacian Gierer-Meinhardt reaction diffusion on a WS network 
with k = 2> and p = 0.1. Concentration of u {top) and v (bottom) on 50 (left) and 500 (right) 
vertex networks. The horizontal line gives the reaction steady state. 



formation followed a Turing mechanism whereby the steady state of the reaction dynamics, 
homogeneous across the vertices, became de-stabilized by the jumps resulting in a non- 
homogeneous pattern across the vertices. The Turing mechanism can also result in pattern 
formation when the network Laplacian is non-symmetric but other pattern formation mech- 
anisms can occur in this case too. Secondly the use of the non-symmetric Laplacian may 
by itself, lead to a build up in concentrations on vertices according to their degree. Thirdly 
if the non-symmetric Laplacian is coupled with reactions that have a finite non-zero steady 
state then the interplay between the two can result in a different pattern across the network. 
The different signatures of these patterns may prove useful in elucidating the underlying 
transport processes in networks when this is not known at the outset. 

The family of Laplacians we have derived should be used when the underlying process may 
be represented by a CTRW. This embodies a large class of diffusive processes on networks. 
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FIG. 10. (Colour Online) Dispersion Relation for Case A {left) and Case B {right) Laplacian 
Gierer-Meinhardt reaction diffusion on a 500 vertex WS network with A; = 3 and p = 0.1. 



• • • • 



CGEiiiiKiimiiiiiiiiminiD 



0.15 



0.10 



U 



10 20 30 40 50 60 

Vertex 



FIG. 11. (Colour Online)Case A Laplacian diffusion on a 63 vertex binary tree network. The 
vertical lines demark sets of vertices with different degrees. 



By considering the underlying stochastic process, the ad-hoc choices of transport operators 
are constrained resulting in a physically consistent model of diffusion on networks. Our 
models of react ion- diffusion on networks are capable of reproducing a wide range of observed 
dynamics, as exemplified by our pattern formation examples. 
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